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ABSTRACT 

We  consider  the  use  of  diagonal  scalings  of  the  Laplacian  matrix  as  precondition- 
ers for  matrices  arising  from  other  second  order  self-adjoint  elliptic  differential  operators. 
It  is  proved  that  if  a  diffusion  operator  with  a  piecewise  constant  but  discontinuous  diffu- 
sion coefficient  is  preconditioned  by  a  diagonal  scaling  of  the  Laplacian,  then,  in  the  limit 
as  the  mesh  size  goes  to  zero,  the  optimal  diagonal  scaling  is  just  the  identity.  This  is  in 
contrast  to  the  case  in  which  the  diffusion  coefficient  is  smoothly  varying,  in  which  case 
numerical  evidence  suggests  that  the  optimal  diagonal  scaling  is  approximately  equal  to 
the  square  root  of  the  diagonal  of  the  matrix. 
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1.  Introduction. 

In  [2]  experiments  were  reported  using  a  numerical  optimization  code  to  determine  the  precondi- 
tioner  of  a  specified  form  which,  for  a  given  coefficient  matrix,  minimized  the  condition  number  of  the 
preconditioned  system.  One  of  the  more  interesting  experiments  involved  finding  the  optimal  diagonal 
scaling  of  the  Laplacian  to  use  as  a  prcconditioner  for  other  second  order  self-adjoint  elliptic  differential 
operators.  Similar  experiments  had  previously  been  carried  out  in  [1],  and  the  use  of  preconditioners  of 
this  form  has  also  been  discussed  in  [6]. 

Let  A;,  be  the  matrix  arising  from  a  finite  element  or  finite  difference  approximation  for  the  problem 

-VaVu=f    inQ  (1.1) 

u=0     ondn, 

where  the  positive  coefficient  a  varies  throughout  the  domain  Q  and  is  bounded  away  from  zero.  Let  A^,  be 
the  Laplacian  matrix  arising  from  the  same  finite  element  or  finite  difference  approximation  for  the  prob- 
lem 

-Au=f    inQ  '  (1.2) 

u=0     onm. 

LetD  be  any  positive  definite  diagonal  matrix.  One  might  consider  using  the  matrix 

M=DA^D  (1.3) 

as  a  prcconditioner  for  the  matrix  Af,  in  an  iterative  algorithm  such  as  the  Chebyshev  or  conjugate  gradient 
method  to  solve  problem  (1.1).  At  each  iteration  it  is  then  necessary  to  solve  a  linear  system  with 
coefficient  matrix  M,  but  such  linear  systems  are  generally  much  easier  to  solve  than  the  original  problem 
with  matrix  A^,.  It  is  trivial  to  invert  the  diagonal  maaix  D,  and,  on  a  uniform  rectangular  grid,  A;,  can  be 
solved  with  a  fast  Poisson  solver.  On  an  irregular  region  A;,  can  be  solved  by  embedding  the  region  in  a 
rectangle  and  using  an  integral  equation  formulation  of  the  problem  [4].  The  number  of  iterations  required 
by  the  Chebyshev  or  conjugate  gradient  algorithms  can  be  bounded  in  terms  of  the  condition  number  of  the 
preconditioned  system,  and  so  one  might  then  ask  what  is  the  best  diagonal  mau-ix  D  to  use  in  (1 .3)  in  order 
to  minimize  this  condition  number.  That  is,  find  a  positive  definite  diagonal  mau-ix  Df,  such  that 
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min  K((DA,D)-M,)  =  k((D,A,D,)-M,)  ,  (1.4) 

D  e  {positive  definite  diagonat  matrices} 


where  K(Af "' A/,)  is  the  ratio  of  the  largest  to  smallest  eigenvalue  of  M~^A^,  or,  the  condition  number  of  the 
symmetrically  preconditioned  matrix,  A/~"^/l;,Af~"^.  This  is  equivalent  to  finding  a  matrix  D^  which 
minimizes 

over  all  positive  definite  diagonal  matrices  D,  since  the  eigenvalues  of  A^'(D~'>\>,D~')  are  the  same  as 
those  of  (DAf,Dy^Af,.  The  problem  was  stated  in  this  second  form  in  [1]. 

In  this  paper  we  prove  a  somewhat  counter-intuitive  result  about  the  optimal  diagonal  scaling  D^, 
when  the  diffusion  coefficient  a  is  piecewise  constant  but  discontinuous.  Both  the  result  and  the  method  of 
proof  became  apparent  from  studying  numerical  results  of  the  optimization  code,  thus  indicating  the  useful- 
ness of  such  a  code  as  a  tool  in  the  study  of  preconditioners.  The  result  is  that  in  the  limit  as  the  mesh  size 
h  goes  to  zero  the  optimal  diagonal  scaling  Df,  approaches  the  identity  (or  a  scalar  multiple  of  the  identity, 
since  scalar  factors  do  not  affect  the  condition  number).  This  is  in  contrast  to  the  case  of  a  smoothly  vary- 
ing diffusion  coefficient  a,  in  which  case  numerical  evidence  suggests  that  the  optimal  diagonal  scaling  D), 
is  approximately  equal  to  the  square  root  of  the  diagonal  of  the  matrix  Af,. 


2.  A  Piecewise  Constant  Diffusion  Coefficient:  Theoretical  Results. 

The  first  theorem  that  we  prove  is  very  general  in  nature,  applying  to  arbitrary  matrices  and  precon- 
ditioners with  a  certain  algebraic  properly.  It  characterizes  a  space  in  which  the  extreme  values  of  the  Ray- 
leigh  quotient  must  be  attained.  The  next  two  theorems  use  this  result  and  apply  to  matrices  arising  from 
specific  forms  of  equation  (1.1),  with  preconditioners  of  the  form  (1.3). 


Theorem  1.  Let  A  and  C  be  two  n  by  n  symmetric  positive  definite  (SPD)  matrices  and  assume  that  certain 
rows  of  C  are  just  scalar  multiples  of  the  corresponding  rows  of  A;  that  is,  there  is  a  nonempty  set  S  such 
that  for  each  ;'  e  S  there  is  a  scalar  c,  such  that 

C.,  =  c,A.j,      Vj  =  \,---,n.  (2.1) 

Then  the  extreme  values  of  the  Raylcigh  quotient  — = —  are  obtained  for  vectors  v  satisfying  either 

V  Cv 

(Av),  =0      Vi£S  (2.2) 

or 

v,=0      VjtS.  (2.3) 


Proof.  Let  w  be  an  arbitrary  vector  and  let  v  be  a  vector  which  satisfies  (2.2)  and  which  matches  w  in  all 
components  outside  of  S.  Such  a  vector  exists  since  A  is  SPD  and  hence  ever)'  principal  submatrix  is  non- 
singular.  The  vector  w  can  be  written  in  the  form 

W  =  V  +  V, 

T    ~      ~r 
where  v  satisfies  (2.3).  Hence  v  Av  =  v  Av  =  0  and  we  have 


w  Aw  =  V  Av  +  V  Av. 
Since  (Cv),  is  also  zero  for  all  i  £  S,  il  also  follows  thai 


w^Cw  =  v^Cv  +  V   Cv. 


Thus  the  Raylcigh  quotienl  is  given  by 


'Aw       vMv  +  V  Av 


w^Cw       v^Cv  +  v  Cv 
and  since  each  term  in  the  numerator  and  denominator  is  nonncgative,  it  satisfies 


v^Av     V  Av         w^Aw  ^  v^Av     v  Av 

<  max 


V  Cv     V  Cv  \       w'Cw  [   v'Cv     V  Cv 

from  which  the  desired  result  follows 


The  spaces  defined  by  (2.2)  and  (2.3)  arc  both  A-orthogonal  and  C-orthogonal  to  one  another,  and  together 
they  span  all  of/?".  The  technique  of  expressing  a  given  vector  as  a  sum  of  vectors  from  each  space  will 
be  used  in  the  proofs  throughout  the  paper. 

Wc  first  use  theorem  1  to  prove  a  result  about  the  one-dimensional  problem 


—-{a—-)=f.      ATE  (0,1) 

ax      ax 
w(0)  =  «(l)  =  0, 


where  the  coefficient  a  (x)  has  the  form 


(2.4) 


a{x)  = 


fli.  ifA:<.5 


a,.  ifx>.5  •      «i.«2>0.    «,^2. 


(2.5) 


Let  A^  be  the  matrix  arising  from  a  continuous  piecewise  linear  finite  element  approximation  for  this  prob- 
lem on  a  uniform  grid  of  size  h.  Assume  that  the  grid  contains  a  node  at  the  point  of  discontinuity  of  a  {x). 
The  matrix  A;,  is  then  given  by 


2a 1  -fli 
-fli 

-a, 
-Ui   2a\     -a  I 

-a I  a]+a2  -a 2 
-a ,     2a n 


A,= 


-0  2  2a 2 


(2.6) 


-4- 


Let  A>,  be  the  one-dimensional  Laplacian  matrix  arising  from  a  continuous  piecewise  linear  finite  element 
approximation  for  the  problem 


uiO)  =  u{l)  =  0, 


(2.7) 


on  the  same  uniform  grid.  The  matrix  A^  is  just  tridi  (-1,2,-1).  We  prove  the  following  theorem. 


Theorem  2.  Let  A^  and  A^  be  as  defined  above  and  let  D^  be  a  positive  definite  diagonal  matrix  which 
satisfies  (1.4).  Then 

(1)     Df,  has  the  form 


d^J, 


do  J, 


where  Ji  ^,  d2,h^  and  df,  are  positive  scalars  and  />,  is  the  identity  of  order  ,  where  h= 

2  n  +  l 

(2)      In  the  limit  as  /!— )0,  these  scalars  approach  each  other;  that  is. 


lim  di  t,  =  lim  ^2  a  =  lim  df,  =  d. 


UD,,  is  any  matrix  of  the  form 


d,,h 


di  J, 


(2.8) 


and  the  positive  scalars  di  >,,  d2,k,  and  dt,  approach  different  limits  as  h-^0  (more  generally,  if  there 
exist  positive   constants   e   and   5   such   that   for  all   h   less   than   6  either     l^,^-^;,l  >  e    or 
\d2,h-dh  I  >  e),  then 

K-((D,A,D,)-'A,)>0(/!-2)      as  h^Q. 


We  prove  this  theorem  through  a  scries  of  lemmas.  For  simplicity  we  drop  the  subscript  h  when  it  is 
clear  which  variables  depend  on  h.  The  point  of  discontinuity  of  a,  x=.5,  is  grid  point  number  — - — ,  which 
we  denote  by  m. 


-5- 


Lemma  2.1.  Let  D  be  any  mairix  of  the  form 


D  = 


d.l 


where  di,  di,  and  d  are  positive  scalars  and  /  is  the  identity  matrix  of  order  m-1.  Define  M  to 
matrix  DAD,  where  A  is  the  Laplacian  matrix.  The  vectors  v  for  which  the  Rayleigh  quotient 
attains  its  extreme  values  satisfy 

(Av),=0,      j  =  l,  •  •  •  ,w-2,  m+2,...,n, 

or,  equivalently, 


n-1-2/ 
n-\-2j 


Vm+1+;  = 


n-1 


7=1,  •  ■  ■  ,m-2 
j  =  l,  ■■■  ,m-2. 


(2.9) 


be  the 
v^Afv 


(2.10) 


(2.11) 


Proof.  Result  (2.10)  follows  from  theorem  1  since  all  rows  of  Af,  except  rows  m-l,  m,  and  w  +  1,  are  just 
scalar  multiples  of  the  corresponding  rows  of  A  and  since  either  of  the  extreme  values,  a\  ld\  or  02/^2. 
which  can  be  taken  on  by  the  Rayleigh  quotient  for  vectors  which  are  zero  outside  fl,...,w-2,w+2,...,«j 
can  also  be  taken  on  for  vectors  satisfying  (2.10). 

By  definition  of  the  matrix  A,  the  equations  (2.10)  are  equivalent  to 


2 

-1 

' 

Vm-2 

Vm-1 

2 

-1 

Vm+2 

-1 

2 

1 

Vm-3 

= 

0 

and 

-1 

2 

-1 

Vm+3 

-1 

-1 

2 

V,  ^ 

0  ^ 

—  1 

1     2 

V„ 

It  is  easy  to  check  that  vectors  of  the  form  (2.11)  are  the  solutions  to  these  equations,  and  so,  the 
equivalence  of  (2.10)  and  (2.1 1). 


Vectors  satisfying  (2.10)  or  (2.1 1)  are  called  discrete  harmonic. 

Lemma_2.2.  (Assertion  (2)  of  the  theorem).  Let  D  be  any  positive  definite  mauix  of  the  form  (2.9).  Ifdi, 
dj,  and  d  approach  different  limits  as  h-^0,  then 


K(M-M)>0(/i~^)     a.v/i-)0, 


where  M=DAD. 


Proof:  For  any  vector  v  satisfying  (2.10)  and  (2.11),  we  can  write 

v^Av  =  v„_i(Av)„_i  +  v„(Av)„  +  v„+,(Av)„+i 

=  v„_iai(2v„_i -v„_i  -  v„)  +  v„((ai+a2)vm-aiv„_i  -a2Vm+i)  + 

n  —  \ 

/T  1  -3  , 

Vm+ia2(2v„+i  -  v„ rv„+i). 

/I  — 1 

After  simplification  this  becomes 


2  2 

v^Av=ai  [(v„-v„_i)^  + -v^_i  ]+a2  [(Vm-Vm+i)^  + rVm+1  ]• 

/i-l  n-\ 


(2.12) 


Similarly,  v  A/v  can  be  written  as 


v^A/v  =  (dv„-diV„.,)^  + r^?^m-i  +('^V„-J2V„+l)^+ r£'2Vm+l. 


n-1 


n-\ 


(2.13) 


Taking  v„_,  =  v„  =  v„+i  =  1  and  dividing  (2.13)  by  (2.12)  gives 
v^A/v 


{d-d,f  +  -^d]  +  {d-dj?  +  -\dl 
n-\  n-\ 


v^Av 


--(ai+a2) 

n-\ 

„-l      m2x{{d,-dfAd2-df} 
2  a  1+02 


>0(h-^). 


Taking  v„_i  =  — ,  v„+i  =  — ,  v„  =  1,  and  dividing  (2.12)  by  (2.13)  gives 
di  di 


r^         ad{\-dld,f+-^{dld,Y]  +  a2[{\-dld2f+^{dld2Y] 
V  Av  n-\  n-\ 


n-\ 


n-\      max  {a^{\ld-\ld^f,  a2{\ld'\ld2f} 


>0{h-^). 


Hence,  by  definition  of  k,  we  have 


K(M-M)  = 


vMv 


v»o  v^A/v 


V  Av 


»o  v^Mv 


vMv 


v*o  v^A/v 


v^A/v 


v»o   v'^^Av 


n-\ 

2 


max^ai(l/rf-l/<ii)^  02(1/^-1/^2)^;      maxf(rfi-^)^  (^2-^)^^  >  ^  ,^-2. 


^1+02 


Lemma  2.3.  (Assertion  (1)  of  ihc  theorem.)  If  D  is  a  matrix  of  the  form  (2.9)  which  satisfies 

min         k((DAD)-U)  =  k((DAD)-'A) 

D  of  the  form  (2.9) 

then  D  also  satisfies 

min  k((DAD)-M)  =  k((DAD)-M). 

Dzipositive  definile  diagonal  matrices} 

Moreover,  any  positive  definite  diagonal  matrix  which  satisfies  this  equation  is  of  the  form  (2.9). 

Proof:  Let  D  =  diag  (5,),  j=l,  •  •  •  ,n  be  any  positive  definite  diagonal  matrix  and  let  Dbe  the  matrix  of  the 
form  (2.9}^  whose  (m-l)",  m'*,  and  {m  +  \y'  diagonal  elements  are  equal  to  those  oi D.  Define  M  =DM> 
and  M  =  DAD.  Let  v  be  a  vector  satisfying  (2.10).  Then  v^Mv  satisfies 

v'''Mv  =  v'''DD'^MD~^Dv  =  w'''Mw, 

where  w=D  £>v  matches  v  in  components  OT-l,m,  and  w  +  1.  AsinTheorem  1,  then,  w  can  be  written  in 
the  form  vv  =  v  +  v,  where  v„_i=v„=v„+i=0,  and  hence  v  Mv  =  v^A/v  =  0.  It  follows  that 

v^Mv  =  w^Mw  =  v^A/v  +  V  Mv  >  v  Mv. 

Since,  by  theorem  1 ,  the  largest  value  of  the  Raylcigh  quotient  -^ —  is  obtained  for  a  vector  v  satisfying 

V  Av 

(2.10),  it  follows  that 

v^Mv  ^         v^Mv  n  T/\\ 

max^= — >max— r — .  (2.14) 

v»o  V  Av        ^*o   V  Av 

Now  let  a  vector  w  be  given  by 

w=D~^Dv, 

where  v  again  satisfies  (2.10)  Then  w^Mw  is  equal  to  v^A/v,  and,  since  the  {m-Vf,  m"',  and  im  +  \y'  ele- 
ments of  vv  match  those  of  v,  we  can  again  write  w  in  the  form  h-  =  v  +  v,  where  v„_]=y„=v„+i=0. 
Hence  v  Av  =  v^Av  =  0  and  we  have 

w^Aw  =  v^Av  +  V  Av  >  v^Av. 

Since,  by  theorem  1,  the  Raylcigh  quotient ^r-  obtains  its  largest  value  for  some  v  satisfying  (2.10),  it 

v^Mv 
follows  that 

max = — >  max ^r— .  (^-i-^) 


8 


From  (2.14)  and  (2.15)  and  the  definition  of  k,  then,  the  desired  result  follows: 

v^Av 


K(Af-'A)<KOW    A)  = 


.0  v^A/v 


max  4^1  <k(M-'a). 


Since  the  inequalities  in  (2.14)  and  (2.15)  are  strict  unless  v  is  zero,  i.e.,  unless  D  is,  itself,  of  the  form 
(2.9),  the  second  part  of  the  lemma  is  also  proved. 


When  the  Laplacian  A^  is  used  as  a  preconditioner  for  A^,,  the  condition  number  of  the  precondi- 
tioned system  is  bounded,  independent  of  h: 


K(A;;Mj<maxf— ,  — ;. 
a-,     a, 


Theorem  2  shows  that,  for  small  h,  the  best  diagonal  scaling  is  close  to  the  identity  (or  a  scalar  multiple  of 
the  identity)  and  so  this  bound  cannot  be  improved  much.  The  wrong  diagonal  scaling  can  greatly  increase 
this  condition  number.  For  larger  values  of  h,  however,  an  appropriate  diagonal  scaling  can  significantly 
reduce  the  condition  number  of  the  preconditioned  system,  as  will  be  demonstrated  in  the  following  sec- 
lion. 

A  similar  result  holds  for  the  two-dimensional  problem 


-{-fia^)+-^{a^))=f,     {x,y)ziO,\)x{0,l) 
ax      ax       ay      ay 

M(x,0)  =  M(x,  \)^uiO,y)  =  u{l,y)^0 


(2.16) 


where  the  coefficient  a  (x,>)  has  the  form 


'^'''^^C^,Z>'s'  ^••^^>«'  '^'^^• 


(2.17) 


Again  let  A^  be  the  matrix  arising  from  a  continuous  piecewise  linear  finite  element  approximation  for  this 
problem  on  a  uniform  triangular  grid  of  size  h,  having  a  mesh  line  at  the  discontinuity,  >=.5.  If  the  natural 
ordering  of  nodes  is  used  then  A,,  has  the  form 


Ah  = 


aj 

-a,I 

-a, I 

-a^I 

-a J  a^T 

-a,I 

-a,l 

-a  2^ 

-a-yl     QnT 


(2.18) 


-9 


where  T=tridii-\A-l)  and  /  is  the  identity  of  order  n  for  a  grid  of  n  by  n  interior  nodes.  Let  A>,  be  the 
two-dimensional  Laplacian  matrix  arising  from  a  continuous  piecewise  linear  finite  element  approximation 
for  the  problem 


-(|-y  +  TT)=/'     (^.>)e(0,1)a:(0,1) 
ax         ay 

u{x,0)  =  u{x,l)  =  u{0,y)  =  u{\,y)  =  0 


(2.19) 


on  the  same  uniform  grid.  The  matrix  A^  is  block  tridi{-IJ,-I).  The  following  theorem  is  proved  very 
similarly  to  the  1-D  case. 

Theorem  3.  Let  A^,  and  A,,  be  as  defined  above  and  let  D>,  be  a  positive  definite  diagonal  matrix  which 
satisfies  (1.4).  Then 

(1)      Dy,  has  the  form 


di.kh 


dj. 


dfhh 


where  di,(,,  ^2,*,  and  d/,  are  positive  scalars,  />,  is  the  identity  of  order  — — — ,  and  /  5  >,  is  the  iden- 


tity of  order  n,  where  h= 


n+\ 


(2)      In  the  limit  as  /i^O,  these  scalars  approach  each  other;  that  is. 


lim  d  1 ;,  =  lim  ^2  A  =  ^"^  ^/i  -  d. 
A-»o      ■         A-»0      '         *-»o 


If  Dl  is  any  matrix  of  the  form 


D.= 


d,,h 


dJ 


diJ) 


(2.20) 


and  the  positive  scalars  d^,,,  1^2./,,  and  df,  approach  different  limits  as  /i— >0  (more  generally,  if  there 
exist  positive   constants  e  and  8  such   that  for  all   h  less  than  8  either     \di_),-dh^  >  e    or 
^di.h-dh  I  >  e),  then 


K{{D,A,D,r'A,)>0{h-^)     as  h^. 
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As  in  the  1-D  case,  we  prove  this  theorem  through  a  series  of  lemmas,  dropping  the  subscript  when  it 
is  clear  which  variables  depend  on  h.  The  matrices  considered  in  the  2-D  case  can  be  thought  of  as  block 

«  4-1 

mauices,  with  n  blocks,  each  of  order  n.  The  subscript  m=  will  denote  the  middle  block,  correspond- 

ing to  the  Hnc  of  discontinuity  in  a.  For  any  n^ -vector  v,  v^  will  denote  the  k"'  block  of  v. 


Lemma  3.1.  Let  D  be  any  matrix  of  the  form 


D  = 


dj 


dl  ■ 


d-,1 


(2.21) 


where  d^,  di,  and  d  are  positive  scalars,  /  is  the  identity  matrix  of  order  — - — ,  and  /  5  is  the  identity 

matrix  of  order  n.  Define  A/  to  be  ihc  matrix  DAD,  where  A  is  the  Laplacian  matrix.  The  vectors  v  for 

which  the  Rayleigh  quotient  -^ —  attains  its  extreme  values  satisfy 
V  Mv 


(Av),=0,      1  =  1,  •••  ,m-2,  OT+2,...,n, 


(2.22) 


Proof:  As  in  the  1  -D  case,  the  result  is  an  immediate  consequence  of  theorem  1 . 


Lemma   3.2.    Let   v   be   a   vector   satisfying   (2.22)   and   let   v„_i    and    v„+i    be   eigenvectors   of 
T  =  rrj'ii/  (-1,4,-1).  Then  each  block  v,  can  be  written  in  the  form 


v,=Y,v„_i,      i  =  \,...,m-2 
v,=Y,v„+i,      i=m+2,...,n 


(2.23) 


for  some  scalars  7,.    If  v„±,    is  the  eigenvector  associated  with  the  smallest  eigenvalue  of  T,  then 
Y.±2=l-0(/i). 


Proof:  Equations  (2.22)  are  equivalent  to 


T  -I 
-I 


-I    7 


"1 

Vm-2 

Vm-1 

1 

V„-3 

= 

0 

and 

I 

.     ^'. 

6 

T  -I 
-/    T 


-1 
-I    7 


Vm+2 
Vm+3 


(2.24) 


We  will  consider  only  the  first  set  of  equations  since  the  second  is  handled  in  exactly  the  same  way. 

Assume  there  is  a  solulion  with  v,  =  Y,Vm-i .  '=1 f^i-l,  for  some  scalars  Y.  Let  ^m-i  correspond  to  an 

eigenvalue  \i  of  T.  Then  equations  (2.24)  become 


H-1 
-1 


lm-2 

c 

-1 

= 

-1     M 

^    Yl  ^ 

.G 

(2.25) 
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This  has  a  unique  solution  for  Yi Ym-2  if  M^2,  and  since  all  eigenvalues  of  T  arc  greater  than  2  this 

condition  holds  and  the  solution  of  (2.24)  is,  indeed,  of  the  form  (2.23). 

Let  Tt  denote  the  tridiagonal  matrix  tridi{-\,\X,-\)  of  order  k  and  let  dcl{'I\)  denote  its  determinant. 
Solving  (2.25)  using  Cramer's  rule  gives 


_  detiT„_^) 
'^"'"^~  dei{T„_2)' 


where  det  (7*)  satisfies 


det{To)=\,     der(r,)  =  n, 
det{T,)  =  ^detiT,.,)-dei{'I\_2\ 


det  a\) 

and  Ti  =  -r— z^ — r  satisfies 

det{Tt.i) 


n  =M. 


ri  =  u. ,     k=2,...,m-2. 

rk-i 


If  carried  out  indefinitely,  this  recurrence  converges  to  a  solution  of  the  equation 

1 


namely, 


2 


and  it  is  easy  to  check  that  after  w -2  =  (9  (/i  ')  steps,  the  ratio  r„_2  is  greater  than  this  limit  by  0(h).  If  |i 
is  the  smallest  eigenvalue  of  T,  then  |i  =  2+0  (/i^),  and  so  r„_2  =  1+0  (h)  and  y„_2  =  1  /r„_2  is  1-0  (/:). 


Lemma  3.3.  (Assertion  (2)  of  the  theorem).  Let  D  be  any  positive  definite  matrix  of  the  form  (2.21).  If 
di,d2,  and  d  approach  different  limits  as  /i— >0,  then 

K{M'^A)>0(.h-^)    as  h^O, 

where  M=DAD. 


Proof:  Let  V  be  a  vector  satisfying  (2.22),  with  v„_i  =  v„  =  v„+i  being  the  eigenvector  of  I,  of  unit  norm, 
corresponding  to  the  smallest  eigenvalue,  |i  =  2+0  (h^).  Then  v^Av  is  given  by 

a  1+6(2 
v^Av  =  a,(|i-I-Y„-2)  +  — z — ^-ai-a2  +  «2(^-l-Ym+2)  =  0(/i) 

while  v^Mv  satisfies 
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v^A/v  =  dl(ii-y„.2}-2did  +  5V  +  dl{[i-y„^2y2dd2  >  {d,-df  +  (d-d^f 
Hence  the  ratio  satisfies 


^>0(/z-'). 
V  Av 


If,  instead  of  having  unit  length,  the  blocks  v„_i  and  v„+i  are  taken  to  have  lengths  — —  and  — — ,  respec- 
lively,  then  we  find 


"1  d-i 

v^A/v  =  d'[0i-Y^-2-l)  +  (M-2)  +  (M-Y».+2-l)]  =  0  {h). 

In  this  case,  then,  we  have 

4^>0(r>). 

and  so  the  condition  number  satisfies 


K{M-^A) 


vMv 


v»o  v'Mv 


v'^Mv 


v»0    vMv 


>0(/j-^). 


Lemma  3.4.  Let  D  be  any  positive  definite  diagonal  matrix  whose  (m-lY',  m'^,  and  (m  +  \y'  diagonal 
blocks  are  just  scalar  multiples  of  the  identity.  Let  D  be  the  matrix  of  the  form  (2.21)  which  matches  D  in 
blocks  m-1,  m,  and  m -1-1.  Then 

K((DAD)-'A)  <  k((5aD)-M). 


Proof:  The  proof  is  analogous  to  that  of  lemma  2.3  in  the  1-D  case. 


Lemma  3.4  shows  that  the  matrix  D  of  the  form  (2.21)  which  minimizes  k((DAD)"'A)  over  all 
matrices  D  of  the  form  (2.21)  also  minimizes  this  quantity  over  all  diagonal  mau-ices  whose  {m-lf,  m**, 
and  {m+iy  diagonal  blocks  are  scalar  multiples  of  the  idcniiiy.  To  show  that  it  minimizxs  this  quantity 
over  all  diagonal  mau^ices,  with  possibly  nonconstant  elements  in  these  blocks,  requires  some  additional 
work.  To  this  end,  wc  prove  the  following  lemma. 


Lemma  3.5.  Let  D  be  any  positive  definite  matrix  of  the  form  (2.21)  and  let  M=DAD.  The  vectors  v  for 

v^Av 
which  the  Raylcigh  quotient  — = —  attains  its  extreme  values  have  blocks  of  the  form 

V  Mv 
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v,  =  a,,v,    1  =  1 n,  (2.26) 

where  Oi ,  •  •  •  ,a„  are  scalars  and  s  is  the  eigenvector  of  T  corresponding  to  the  smallest  eigenvalue. 


Proof:  Let  5  be  the  matrix  whose  columns  are  the  eigenvectors  of  7",  and  let  6  be  the  diagonal  matrix  of 
eigenvalues,  so  that  we  have  TS  =  56,  S^S  =  SS^  =  I.  Define  [/  to  be  the  block  diagonal  matrix  whose 
diagonal  blocks  are  all  equal  to  S.  Then  U^AU  is  of  the  form  (2.18)  with  T  replaced  by  8  and  U^MU  is 
D*  (block  tridi(-I,Q-I))*D.  Let  P  be  a  permutation  matrix  such  that  the)""  column  of  the  i'*  block  ofa 
matrix  B  is  the  /'"■  column  of  the  j'''  block  of  the  mau-ix  BP.  Then  multiplying  by  P^  on  the  left  and  by  P 
on  the  right,  the  above  matrices  become  block  diagonal  with  iridiagonal  blocks  given  by 


and 


ajB,  -a  I 
-ai  -a, 

-a  1   a  1 0,        -a  I 
ai+a2 


-^2 

-02  -02 

-0  2    0  2^, 


-d'  -d} 


-d\   d\Q,    -d^d 


-dd2   ^2^1   -c/j 


-di 


-dl  dl^^ 


respectively.  Clearly,  the  exu-eme  values  of  the  ratio  w^P^U^AUPw/w^P^U^MUPw  are  attained  for  vec- 
tors H'  with  a  single  nonzero  block,  corresponding  to  the  smallest  eigenvalue  9,.  The  vector  v  =  UPw,  then, 
is  an  extreme  vector  for  the  Rayleigh  quotient  v^Av/v^A/v,  and  each  block  is  a  scalar  multiple  of  the  eigen- 
vector 5  of  r corresponding  to  the  smallest  eigenvalue. 


Lemma  3.6.  (Assertion  (1)  of  the  theorem.)  If  D  is  a  matrix  of  the  form  (2.21)  which  satisfies 


min         K((DAD)-'A)=K((DAD)-U) 

D  of  the  form  (Z2\) 


(2.27) 


then  D  also  satisfies 


min  k((DAD)-M)  =  k((Z)AD)-M). 

Dzfpositive  defirute  diagonal  matrices} 


(2.28) 
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Moreover,  any  positive  definite  diagonal  matrix  which  satisfies  (2.28)  is  of  the  form  (2.21). 


Proof:  Let  D  =  diag{Di,  ■  ■  ■  ,D„)  be  any  positive  definite  diagonal  matrix.  Let  D  be  the  matrix  of  the 
form  (2.21)  whose  {m-\y',  m"',  and  im  +  \y'  block  coefficients  are 

di=s'^D„_iS,     d  =  s'^D„s,     d2=s'''D„^_iS, 

where  5  is  the  eigenvector  of  T  corresponding  to  the  smallest  eigenvalue  p..    Define  M=DAD  and 
M  =  DAD.  Let  v  be  a  vector  satisfying  (2.22)  and  (2.26).  Then  v^Afv  satisfies 

v^Mv  =  v'^DD'^MD'^Dv  =  w'^Mw, 

where  w=D    Dv.  The  vector  w  can  be  written  in  the  form  v+v,  where  v  =  (D    D-I)v.  Because  of  the 
choice  of  di,  di,  and  d,  we  have 

V  Mv  =  a„_i5^(^'D„_i  -I)^{d^a„-i\is  -dia„_2X  -dida„s)  + 
a„s  {d    D„-I){d  a„\]LS- dida„-is -dd20.„+xs)  + 
a„+i5^(rfi'£»„+i  - 7)^(^2 a„+i|jj  -d2a„+2.s  -dd20.„s) 
=  0. 

It  follows  that 

v^A/v  =  w^Mw  =  v^A/v  +  V  Mv  >  v^A/v. 


v^Mv 
Since,  by  theorem  1  and  lemma  3.5,  the  largest  value  of  the  Rayleigh  quotient  -^ —  is  obtained  for  a  vec- 

vMv 

lor  V  satisfying  (2.22)  and  (2.26),  it  follows  that 

v^Mv  ^         v^Mv  ,,  ^„, 

max^^ — >max^^ — .  (2.29) 

v»o   V  Av       v*o   V  Av 

Now  let  a  vector  w  be  given  by 


where  v  again  satisfies  (2.22J  and  (2.26).  Then  w^Mw  is  equal  to  v^A/v,  and  again  we  can  write  >v  in  the 
form  w  =  V  +  V,  where  v  =  (D    D-I)v.  We  now  have 

V  Av  =  a„_i.?  (rfiD„_i-/)'ai(a„_i|ij  -a„_25  -a„s)  + 

_-_i  a\+a2 

a„s'(dD„-Iy{ — - — a„\is  -  a^a„^]S  -  ajCL^+^s)  + 


am+li^('^2^m  +  l-/)^a2(am  +  lU^  -  O-mS  -  o.„^2s). 
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which  can  be  written  in  the  form 

v^Av  =  v^_i  (C/lCv)„_,  +  vl{CACv)„  +  v^i  (CACv)„^i , 

where  C  is  a  diagonal  matrix  whose  diagonal  elements  are  one  except  in  blocks  m-\,  m,  and  m  +  \,  where 
they  are 

respectively.  Because  of  the  choice  of  d],  d2,and  d,  we  know  that  the  quantities  under  the  square  roots  are 
nonnegalive,  since 

;  =  1  /  =  1  j=\  j=\  t  =  i  d,,,       d,j 

>  =  1  /  =  1  t=l  ;=1 

Hence  v  Av  is  nonnegative  and  so  we  have 

w^Aw>v^Av+v  Av>v^Av. 

Since,  by  theorem  1  and  lemma  3.5,  the  Raylcigh  quotient r—  obtains  its  largest  value  for  some  v  satis- 

fying  (2.22)  and  (2.26),  it  follows  that 

w'^Aw               v^Av  ,-  ,„, 

max = — >  max ^.  (2.jU) 

Combining  (2.29)  and  (2.30),  we  obtain  the  desired  result: 

K(A/"' A)  <  K(Af  "' A)  <  k(M''a). 

Since  the  inequalities  in  (2.29)  and  (2.30)  are  strict  unless  v  is  zero,  i.e.,  unless  D  is,  itself,  of  the  form 
(2.21),  the  second  part  of  the  lemma  is  also  proved. 


While  our  primary  interest  has  been  in  diagonal  scalings  of  the  Laplacian,  it  should  be  noted  that  the 
proofs  of  lemmas  2.3  and  3.6  make  no  use  of  the  assumption  that  D  is  diagonal  outside  of  positions 
(blocks)  m-1,  m,  and  m+l.  They  can  therefore  be  generalized  to  the  following  result: 

Corollary.  For  the  1-D  problem,  the  matrix  D^  of  theorem  2  minimizes  k((£[A;,£;,)"' A  J  over  all  matrices 
Eh  who.se  three  center  rows  and  columns  (w-1,  m,  and  m  +  \)  have  nonzeros  only  on  the  diagonal.  For  the 
2-D  problem,  the  mau-ix  D^,  of  theorem  3  minimizes  k((£'[A;,£J"'A^)  over  all  matrices  Ef,  whose  3;j  center 
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rows  and  columns  (n{m-\)  through  n(m  +  \))  have  nonzeros  only  on  the  diagonal  and  these  nonzeros  are 
positive.  More  generally,  it  minimizes  this  quantity  over  all  matrices  whose  3n  center  rows  and  columns 
have  nonzeros  only  in  the  n  by  n  diagonal  blocks  and  for  which  these  diagonal  blocks,  £„_i ,  E„,  and  £„+i , 
have  the  property  that 

(s'^E,s)is'^E-^s)>\,     i=m-\,m,m  +  \, 

where  s  is  the  eigenvector  of  T  corresponding  to  the  smallest  eigenvalue. 


3.  Numerical  Results. 

For  a  given  matrix  Af,  an  optimization  code  can  be  used  to  determine  numerically  the  optimal 
preconditioner  of  the  form  (1.3).  A  particularly  efficient  technique  for  solving  this  type  of  optimization 
problem  was  developed  by  Overton  [5].  Experiments  with  this  code  were  reported  in  [2].  The  code  uses  a 
variant  of  Newton's  method  to  determine  the  matrix  A/  of  a  specified  form  (e.g.,  form  (1.3))  for  which  the 
spectral  radius 

p(/-A/-M,) 

is  minimal.  It  was  shown  in  [2]  that  this  same  matrix  M  (or  any  scalar  multiple  of  M)  also  minimizes  the 
condition  number  k(A/''A>,),  provided  the  set  over  which  the  minimization  is  being  performed  contains  all 
positive  scalar  multiples  of  its  members,  which  it  does  in  this  case. 

In  the  following  experiment,  the  matrix  Af,  was  taken  to  be  the  matrix  arising  from  a  continuous 
piecewise  linear  finite  element  approximation  on  a  uniform  grid  of  size  h  for  the  one  dimensional  problem 
(2.4)-(2.5),  where 

ai  =  l,     02  =  100-  0-1) 

The  optimization  code  was  run  to  determine  the  optimal  preconditioner  of  the  form  (1.3).  The  diagonal 
matrix  D^  determined  by  the  code  was  always  of  the  form  (2.9),  as  Theorem  2  shows  it  must  be.  (In  fact, 
this  observation  of  the  numerical  results  led  to  the  statement  and  proof  of  Theorem  2!)  The  values  of  the 
diagonal  elements  d^,  d2,  and  d  and  the  condition  number  k  of  the  optimally  preconditioned  system  are 
listed  in  Table  1  for  various  grid  sizes.  The  actual  matrix  D;,  returned  by  the  code  has  been  multiplied  by  a 
scalar  so  that  d  is  1 . 

The  largest  problem  size  that  the  optimization  code  was  able  to  handle  at  the  time  of  these  tests  was 
about  n=:225,  but  it  is  currently  being  modified  to  handle  larger  matrices.  For  this  value  of  n,  the  element 
dx  and  the  condition  number  k  have  reached  only  about  half  of  their  asymptotic  limit.  It  is  not  at  all  clear 
from  the  numerical  results  alone  that  there  is  an  asymptotic  limit,  since  this  is  approached  only  for  much 
smaller  values  of  h.  This  leads  one  to  question  the  relevance  of  asymptotic  results  such  as  that  in  Theorem 
2,  since  for  typical  size  problems  they  may  not  be  approached.  It  is  interesting  to  note  that  for  all  problem 
sizes  the  scalar  d  is  equal  to  di,  the  diagonal  element  corresponding  to  the  subregion  with  the  larger  diffu- 
sion coefficient. 
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1//I 

d, 

d2 

d 

K{M-'A,) 

10 

.191 

1.000 

1.22 

26 

.272 

1.000 

14.69 

50 

.340 

1.000 

22.89 

82 

.396 

1.000 

31.10 

122 

.444 

1.000 

39.08 

226 

.542 

1.000 

53.50 

Tabic  1.  Scalars  Defining  the  Optimal  Preconditioner  and  Condition  Number  of  the  Preconditioned  System 
for  the  One  Dimensional  Problem  (3.1). 


Table  2  shows  numerical  results  for  a  similar  two  dimensional  problem: 


-(^-a^-  + 3-03-)=/    w  (0,l)x(0.1), 
dx     ax       ay     ay 


u{x,0)  =  uix,l)  =  u  {Q,y)  =  u  {\,y)  =  0, 


where 


a{x,y)- 


1,       )'<.5 
100,     >'>.5 


(3.2) 


The  matrix  A^  was  again  derived  from  a  continuous  piecewise  linear  finite  element  approximation  on  a  uni- 
form grid  of  size  h.  The  optimization  code  was  used  to  find  the  diagonal  matrix  D^  for  which 

K((D,A,D,)-U,) 

is  minimal,  where  A/,  is  the  5-point  Laplacian. 

The  mauix  Z)^  was  again  observed  to  have  the  form 


D,= 


dJ 


dl 


dJ 


where  d^,  di,  and  d  are  scalars,  /  is  the  identity  corresponding  to  the  subregjon  (0,1)j:(0,.5)  or 
(0,l)x(.5,l),  and/  5  is  the  identity  on  the  dividing  line, >>=. 5.  The  values  di,  ^2.  and  das  well  as  the  con- 
dition number  k  of  the  optimally  preconditioned  system  are  hsled  in  Table  2  for  various  grid  sizes.  Again, 
since  the  optimization  code  was  not  able  to  handle  very  large  matrices,  the  asymptotic  behavior  of  the  sys- 
tem cannot  be  deduced  from  the  numerical  results  alone,  but  is  established  by  theorem  3. 
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1/h 

d^ 

d2 

d 

K{M-'A,) 

4 

.146 

1.251 

2.30 

6 

.152 

1.143 

3.23 

8 

.158 

1.076 

4.07 

10 

.164 

1.034 

4.85 

16 

.187 

1.000 

6.93 

Table  2.  Scalars  Defining  the  Optimal  Preconditioner  and  Condition  Number  of  the  Preconditioned  System 
for  the  Two  Dimensional  Problem  (3.2). 


In  contrast  to  the  above  results,  Table  3  shows  results  for  the  problem  (2.1)  where  a  (x)  is  given  by 

a{x)  =  m+x'^.  (3.3) 

Although  the  total  variation  in  a{x)  over  the  interval  (0,1)  is  approximately  the  same  as  that  in  (3.1) 
(^max/^mm  =  101),  it  now  varics  smoothly.  The  diagonal  matrix  D>,  returned  by  the  optimization  code  is 
now  very  nearly  equal  to  the  square  root  of  the  diagonal  of  >4>,.  Table  3  shows  the  largest  and  smallest  ratio 
between  the  square  of  a  diagonal  element  of  D^  and  the  corresponding  element  oi  A^,.  D^  has  been  multi- 
plied by  a  scalar  so  that  its  center  element  is  equal  to  the  square  root  of  the  corresponding  diagonal  element 
oi  Ah.  In  this  case,  then,  the  optimal  diagonal  matrix  D^  is  not  of  the  form  (2.9)  and  it  does  not  appear  to 
approach  the  identity  as  /z— >0.  Rather  it  appears  to  approach  the  square  root  of  the  diagonal  of  A^. 


llh 

max  oil  A,, 

min  dIiA,, 

K(A/-'A,) 

10 

1.01 

.98 

1.17 

26 

1.00 

.98 

1.26 

50 

1.00 

.99 

1.28 

Table  3.    Ratios  of  Diagonal  Elements  for  the  Optimal  Preconditioner  and  Condition  Number  of  the 
Preconditioned  System  for  the  Problem  (3.3). 


4.  Further  Discussion. 

In  the  case  of  a  smoothly  varying  diffusion  coefficient  a,  the  differential  operator  V  •  (aVw)  can  be 
written  in  the  form 


aAu  +  Va-Vu 


(4.1) 


Consider  the  equation  Au  =-/.   Making  a  change  of  variable  from  u  lo  v  =  a  "^«  and  multiplying  this 
equation  by  a"^  gives 


a"2A(a"2v)  =  aAv  +  Va-Vv+a^'^{Aa^'^)v  =  -a^'^f 


(4.2) 


The  matrix  M  =  DAt,D,  where  D  is  the  square  root  of  ihe  diagonal  of  A^,  represents  the  differential  opera- 
tor in  (4.2),  with  the  same  homogeneous  Dirichlcl  boundary  conditions  as  the  original  problem.  Since  this 
is  a  second  order  sclf-ad joint  operator,  it  follows  that  using  A/  as  a  preconditioner  for  A^  gives  a  condition 
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number  for  the  preconditioned  system  that  is  0(1),  independent  of  the  mesh  size  [3].  Since  the  leading 
terms  of  the  differential  operator  in  (4.2)  match  those  in  (4.1),  it  is  perhaps  not  surprising  that  this  is  a 
near-optimal  diagonal  scaling. 

When  the  coefficient  a  is  discontinuous  or  continuous  but  not  differentiablc,  there  is  no  such  analogy 
between  the  preconditioner  and  a  differential  operator  whose  leading  terTn(s)  match  those  of  the  original 
equation.  In  this  case,  a  discontinuous  diagonal  scaling  of  the  Laplacian  does  not  represent  a  second-order 
self-adjoint  elliptic  operator  and,  as  theorems  2  and  3  show  for  a  specific  problem  class,  the  condition 
number  of  the  matrix  preconditioned  in  this  way  may  become  infinite  as  h  goes  to  zero. 
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